## Script to run model and plots of rebel speech by rebellion
## Replicates Table 4 and Figure 3 in paper
## JBS: final revision revised 3 August 2017


rm(list=ls())

library(plyr)
library(lme4)
library(stargazer)
library(arm)
library(mvtnorm)

##### SET YOUR PATH HERE ##############
##########################################

path <- "~/Dropbox (Personal)/Rebel Summaries/APSR_SKLLO_Repfiles/"

setwd(path)

load("CalcData/rebeldta0510.Rda")
load("CalcData/rebeldta9297.Rda")


mod1 <- glmer(matrix(c(Rebspeech , totspeech), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.9297[all.dta.9297$Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod2 <- glmer(matrix(c(Rebspeech , totspeech), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.9297[all.dta.9297$Party=="Con",],control=glmerControl(optimizer="bobyqa"))

mod3 <- glmer(matrix(c(Rebspeech , totspeech), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.0510[all.dta.0510$Party=="Lab",],control=glmerControl(optimizer="bobyqa"))

mod4 <- glmer(matrix(c(Rebspeech , totspeech), ncol=2) ~ gov*q75+Majority+leader+start_year+(1|PubWhipID), family=binomial(link=logit), data= all.dta.0510[all.dta.0510$Party=="Con",],control=glmerControl(optimizer="bobyqa"))

summary(mod1)
summary(mod2)
summary(mod3)
summary(mod4)


stargazer(mod1,mod2,mod3,mod4,type="latex",title="Effect of Government Status on Rebellious Speech: Fractional Logit Models with Random Effects", align=T, column.labels=c("Labour","Conservative","Labour","Conservative"),dep.var.labels=c("1992--2001","2005--2010"),style="ajps")


#### Create Plots of Simulations of Proportion of Rebel Speeches
#### Simulate backbenchers with mean victory margin (majority)
# Set number of simulations
sims<-1000

# Take random draws
dist<-rmvnorm(sims, fixef(mod1), as.matrix(vcov(mod1)))  # LABS
dist2<-rmvnorm(sims, fixef(mod2), as.matrix(vcov(mod2))) # CONS
dist3<-rmvnorm(sims, fixef(mod3), as.matrix(vcov(mod3))) # LABS
dist4<-rmvnorm(sims, fixef(mod4), as.matrix(vcov(mod4))) # CONS

labdta9297 <- all.dta.9297[all.dta.9297$Party=="Lab",]
condta9297 <- all.dta.9297[all.dta.9297$Party=="Con",]
labdta0510 <- all.dta.0510[all.dta.0510$Party=="Lab",]
condta0510 <- all.dta.0510[all.dta.0510$Party=="Con",]

# Z is indicator of in gov or not
z<-c(0,1)

calc<-length(z)

lo1<-numeric(calc)
hi1<-numeric(calc)
pe1<-numeric(calc)
lo2<-numeric(calc)
hi2<-numeric(calc)
pe2<-numeric(calc)
lo3<-numeric(calc)
hi3<-numeric(calc)
pe3<-numeric(calc)
lo4<-numeric(calc)
hi4<-numeric(calc)
pe4<-numeric(calc)
lo5<-numeric(calc)
hi5<-numeric(calc)
pe5<-numeric(calc)
lo6<-numeric(calc)
hi6<-numeric(calc)
pe6<-numeric(calc)
lo7<-numeric(calc)
hi7<-numeric(calc)
pe7<-numeric(calc)
lo8<-numeric(calc)
hi8<-numeric(calc)
pe8<-numeric(calc)

for(i in 1:calc){
  pp<-numeric(sims)
  pp2<-numeric(sims)
  pp3<-numeric(sims)
  pp4<-numeric(sims)
  pp5<-numeric(sims)
  pp6<-numeric(sims)
  pp7<-numeric(sims)
  pp8<-numeric(sims)
  
  for(j in 1:sims){
    #Labour (92-97); Moderate
    pp[j] <-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*0+dist[j,7]*0*z[i]+dist[j,4]*mean(labdta9297$Majority, na.rm=T)+dist[j,6]*mean(labdta9297$start_year, na.rm=T)
    #Labour (92-97); extreme
    pp2[j]<-dist[j,1]+dist[j,2]*z[i]+dist[j,3]*1+dist[j,7]*1*z[i]+dist[j,4]*mean(labdta9297$Majority, na.rm=T)+dist[j,6]*mean(labdta9297$start_year, na.rm=T)
    #Conservative (92-97); Moderate
    pp3[j]<-dist2[j,1]+dist2[j,2]*z[i]+dist2[j,3]*0+dist2[j,7]*0*z[i]+dist2[j,4]*mean(condta9297$Majority, na.rm=T)+dist2[j,6]*mean(condta9297$start_year, na.rm=T)
    #Conservative (92-97); Extreme
    pp4[j]<-dist2[j,1]+dist2[j,2]*z[i]+dist2[j,3]*1+dist2[j,7]*1*z[i]+dist2[j,4]*mean(condta9297$Majority, na.rm=T)+dist2[j,6]*mean(condta9297$start_year, na.rm=T)
    #Labour (05-10); Moderate
    pp5[j]<-dist3[j,1]+dist3[j,2]*z[i]+dist3[j,3]*0+dist3[j,7]*0*z[i]+dist3[j,4]*mean(labdta0510$Majority, na.rm=T)+dist3[j,6]*mean(labdta0510$start_year, na.rm=T)
    #Labour (05-10); extreme
    pp6[j]<-dist3[j,1]+dist3[j,2]*z[i]+dist3[j,3]*1+dist3[j,7]*1*z[i]+dist3[j,4]*mean(labdta0510$Majority, na.rm=T)+dist3[j,6]*mean(labdta0510$start_year, na.rm=T)
    #Conservative (05-10); Moderate
    pp7[j]<-dist4[j,1]+dist4[j,2]*z[i]+dist4[j,3]*0+dist4[j,7]*0*z[i]+dist4[j,4]*mean(condta0510$Majority, na.rm=T)+dist4[j,6]*mean(condta0510$start_year, na.rm=T)
    #Conservative (05-10); Extreme
    pp8[j]<-dist4[j,1]+dist4[j,2]*z[i]+dist4[j,3]*1+dist4[j,7]*1*z[i]+dist4[j,4]*mean(condta0510$Majority, na.rm=T)+dist4[j,6]*mean(condta0510$start_year, na.rm=T)
  }
  
  pe1[i]<-invlogit(mean(pp))
  lo1[i]<-invlogit(quantile(pp, 0.025))
  hi1[i]<-invlogit(quantile(pp, 0.975))
  
  pe2[i]<-invlogit(mean(pp2))
  lo2[i]<-invlogit(quantile(pp2, 0.025))
  hi2[i]<-invlogit(quantile(pp2, 0.975))
  
  pe3[i]<-invlogit(mean(pp3))
  lo3[i]<-invlogit(quantile(pp3, 0.025))
  hi3[i]<-invlogit(quantile(pp3, 0.975))
  
  pe4[i]<-invlogit(mean(pp4))
  lo4[i]<-invlogit(quantile(pp4, 0.025))
  hi4[i]<-invlogit(quantile(pp4, 0.975))
  
  pe5[i]<-invlogit(mean(pp5))
  lo5[i]<-invlogit(quantile(pp5, 0.025))
  hi5[i]<-invlogit(quantile(pp5, 0.975))
  
  pe6[i]<-invlogit(mean(pp6))
  lo6[i]<-invlogit(quantile(pp6, 0.025))
  hi6[i]<-invlogit(quantile(pp6, 0.975))
  
  pe7[i]<-invlogit(mean(pp7))
  lo7[i]<-invlogit(quantile(pp7, 0.025))
  hi7[i]<-invlogit(quantile(pp7, 0.975))
  
  pe8[i]<-invlogit(mean(pp8))
  lo8[i]<-invlogit(quantile(pp8, 0.025))
  hi8[i]<-invlogit(quantile(pp8, 0.975))
  
  print(i)
}


#92-2001
pdf("FinalPlots/TotRebSpeech92_2001.pdf")
#Labor Loyalists
plot(c(0.15, 0.65), pe1, pch=15, xlim=c(0, 1), cex=1.5, ylim=c(0, 0.2), ylab="Predicted Number of Speeches", xaxt="n", col="grey80", xlab="")
segments(0.15, lo1[1], 0.15, hi1[1], lwd=2, col="grey80")
segments(0.65, lo1[2], 0.65, hi1[2], lwd=2, col="grey80")
grid()

#Labor Rebels
points(c(0.20, 0.70), pe2, pch=16, cex=1.5, col="grey80")
segments(0.20, lo2[1], 0.20, hi2[1], lwd=2, col="grey80")
segments(0.70, lo2[2], 0.70, hi2[2], lwd=2, col="grey80")

#Conservative Loyalists
points(c(0.25, 0.75), pe3, pch=15, cex=1.5, col="grey20")
segments(0.25, lo3[1], 0.25, hi3[1], lwd=2, col="grey20")
segments(0.75, lo3[2], 0.75, hi3[2], lwd=2, col="grey20")

#Conservative Rebels
points(c(0.30, 0.80), pe4, pch=16, cex=1.75, col="grey20")
segments(0.30, lo4[1], 0.30, hi4[1], lwd=2, col="grey20")
segments(0.80, lo4[2], 0.80, hi4[2], lwd=2, col="grey20")

axis(1, at=c(0.225, 0.725), labels=c("Opposition", "Government"), tick=TRUE)
legend("topleft", inset=0.05, c("Labour", "Conservative", "Loyalists", "Rebels", "95% CI"), cex=0.8,  pch=c(NA, NA, 15, 16, NA), lty=c(1,1,NA,NA,1), col=c("grey80", "grey20", "black", "black", "black"))
dev.off()

#2005-2010
pdf("FinalPlots/TotRebSpeech05_2010.pdf")
#Labor Loyalists
plot(c(0.15, 0.65), pe5, pch=15, xlim=c(0, 1), cex=1.5, ylim=c(0, 0.2), ylab="Predicted Number of Speeches", xaxt="n", col="grey80", main="Proportion of Rebel Speeches\n by Frequency of Rebellion and Party Status", xlab="")
segments(0.15, lo5[1], 0.15, hi5[1], lwd=2, col="grey80")
segments(0.65, lo5[2], 0.65, hi5[2], lwd=2, col="grey80")
grid()

#Labor Rebels
points(c(0.20, 0.70), pe6, pch=16, cex=1.5, col="grey80")
segments(0.20, lo6[1], 0.20, hi6[1], lwd=2, col="grey80")
segments(0.70, lo6[2], 0.70, hi6[2], lwd=2, col="grey80")

#Conservative Loyalists
points(c(0.25, 0.75), pe7, pch=15, cex=1.5, col="grey20")
segments(0.25, lo7[1], 0.25, hi7[1], lwd=2, col="grey20")
segments(0.75, lo7[2], 0.75, hi7[2], lwd=2, col="grey20")

#Conservative Rebels
points(c(0.30, 0.80), pe8, pch=16, cex=1.75, col="grey20")
segments(0.30, lo8[1], 0.30, hi8[1], lwd=2, col="grey20")
segments(0.80, lo8[2], 0.80, hi8[2], lwd=2, col="grey20")

axis(1, at=c(0.225, 0.725), labels=c("Opposition", "Government"), tick=TRUE)
legend("topleft", inset=0.05, c("Labour", "Conservative", "Loyalists", "Rebels", "95% CI"), cex=0.8,  pch=c(NA, NA, 15, 16, NA), lty=c(1,1,NA,NA,1), col=c("grey80", "grey20", "black", "black", "black"))
dev.off()
